Writing functions in R

Palaeoverse workshop 2023

Function basics

What is a function

takes input –> does something –> returns output

mean(c(1, 2, 3))
[1] 2


A function needs a name, arguments in (), and a body in {} 

subtract <- function(arg1, arg2) { 
  arg1 - arg2 
} 



subtract(2, 1)
[1] 1

Why do we need functions

  • Readability
  • Organisation
  • Modularity
  • Reusability

Imagine calculating the mean without standard functions like mean or sum:

  data <- c(1,2,3)
  total <- 0
  count <- 0
  for (value in data) {
    total <- total + value
    count <- count + 1
  }
  total/count
[1] 2

Arguments

Arguments need to be provided in the correct order, or specified by name: 

subtract(2, 1)
[1] 1


subtract(1, 2)
[1] -1


subtract(arg2 = 1, arg1 = 2)
[1] 1

Default values

Make function use more convenient, can hide complexities.

subtract(2)
Error in subtract(2): argument "arg2" is missing, with no default


This will work if we set a default for arg2:

subtract <- function(arg1, arg2 = 1) {
  arg1 - arg2
}


subtract(2)
[1] 1


Ellipsis (‘…’)

Additional, optional arguments can be allowed by using ‘…’ as the last argument:

my_plot <- function(arg1, arg2, ...) {
  plot(arg1, arg2, ...)
}


my_plot(2, 1, col = "red", pch = 17, cex = 2)

Return

A function generally should return something, but this does not:

subtract <- function(arg1, arg2) {
  result <- arg1 - arg2
}
subtract(2,1)


Return explicitly with return, or place return value at the end of the function:

subtract <- function(arg1, arg2) {
  result <- arg1 - arg2
  return(result)
}
subtract(2,1)
[1] 1


Return multiple objects

just_return <- function(arg1, arg2) {
return(arg1)
return(arg2)
}
just_return(2, 1)
[1] 2

This did not work as intended. R functions only return one object. Instead use lists or other data structures:

just_return <- function(arg1, arg2) {
return(c(arg1, arg2))
}
just_return(2, 1)
[1] 2 1

Binary operators

Standard function syntax:

sum(c(1,2))
[1] 3

Operator syntax:

1 + 2
[1] 3

Most binary operators come in %:

3 %in% c(1,2,3)
[1] TRUE


Custom binary operators – let’s define an operator for “not in”:

`%!in%` <- function(x, y) !(x %in% y)
3 %!in% c(1,2,3)
[1] FALSE

Control structures – if

if a condition is true, do something.

if (1 + 1 == 2) print("True")
[1] "True"


add_or_subtract <- function(arg1, arg2, operation) {
 if (operation == "add") {
   result <- arg1 + arg2
 }
 if (operation == "subtract") {
   result <- arg1 + arg2
 }
 result
}
add_or_subtract(2,1,"add")
[1] 3

Control structures – else

else instructs what to do when the if condition is not met.

if (1 + 1 == 3) print("True") else print("False")
[1] "False"


add_or_subtract <- function(arg1, arg2, operation) {
 if (operation == "add") {
   result <- arg1 + arg2
 } else {
   result <- arg1 - arg2
 }
 return(result)
}
add_or_subtract(2,1,"subtract")
[1] 1

Control structures – switch

Instead of many if and else statements, try switch

fossil_description <- function(fossil) {
 switch(fossil,
  ammonite = "coiled shell",
  Tyrannosaurus = "serrated teeth",
  Lepidodendron = "scaly bark", 
  "not a fossil"
 )
}
fossil_description("Tyrannosaurus")
[1] "serrated teeth"
fossil_description("Lewis")
[1] "not a fossil"

Control structures – for loops

Loops are used for repeating similar actions multiple times. for loops iterate over a set of values. The iterator (i) changes with every iteration of the loop:

for(i in c(1,2,3)) print(i)
[1] 1
[1] 2
[1] 3

To generate sequences of integers, we can use seq_len. Let’s make a function:

print_repetitions <- function(n) {
 for (i in seq_len(n)) { 
   print(i)
 }
}
print_repetitions(2)
[1] 1
[1] 2

Control structures – while loops

while loops repeat a task until a condition is no longer met.

add_until_4 <- function(x) {
  while(x < 4) {
    x <- x + 1
    print(x)
  }
}
add_until_4(1)
[1] 2
[1] 3
[1] 4

Exercise 1 - Function for latitudinal binning

Create a function that can sort a data.frame into latitudinal bins. That is, we want a new column that identifies the bin of each entry of the data set. As an exemplary data set, we can use the reefs data from palaeoverse.

If you are new to writing R functions, try a simpler function that can sort data into the northern and southern hemisphere.

Here is what the result may look like when sorted into hemispheres:

reefs[72:73, c("name", "lat", "hemisphere")]
          name     lat hemisphere
72 Begunjscica 46.4333      north
73     W-Ceram -3.2500      south

Exercise 1 - My solution

Best practices

Naming style

Give variables and functions consistent names. These are the two most common styles:

snake_case

used by Tidyverse’s style guide

CamelCase

used by Google’s R style guide

Clean coding

Common practices:

  • Avoid renaming existing functions and variables
mean <- mean(c(1, 2, 3))
mean
[1] 2
mean(c(1, 2, 3))
[1] 2

R is clever, in this case this still works.

  • Use <- for assignment, not =
data = c(1, 2, 3)
data <- c(1, 2, 3)

Example of a detailed style guide: Tidyverse’s style guide

Comments

  • document purpose and usage of code
  • explain complex / non-intuitive code
# Calculate convex hull
tmp <- tmp[chull(x = tmp[, lng], 
                 y = tmp[, lat]), c(lng, lat)]
  • organise code into sections
#=== Set-up ===
  unique_taxa <- unique(occdf[, name])
  # Order taxa
  unique_taxa <- unique_taxa[order(unique_taxa)]
#=== convex hull  ===

Good code shouldn’t need a large amount of comments - but comment enough that you can still use your code two years later.

Comments

Add some general information in the beginning of a large R script.

### Change point regression analysis
### July 2021
### Kilian Eichenseer
###
### Bayesian algorithm for finding a change point in 
### the linear relationship between two variables. 
### Uses JAGS (https://mcmc-jags.sourceforge.io/).

### Generate data
set.seed(10) 
n <- 60 # total number of data points

Even better: add formal documentation.

Documentation

?mean

roxygen2

  • R package, install with install.packages("roxygen2")
  • Used to create documentation for R packages (functions, data, …)
  • Start every line of documentation with ’#

This could generate documentation for the subtract function from earlier:

#' Subtraction
#' 
#' Subtracts `arg2` from `arg1`
#'
#' @param arg1 `Numeric`. First argument.
#' @param arg2 `Numeric`. Second argument.
#' @return A `numeric` containing the difference 
#' between `arg1` and `arg2`.
#' @examples
#' subtract(2,1)

?subtract

Exercise 2 - Document your function

  • Make sure your code follows the same style throughout
  • Add in-line comments if necessary
  • Create documentation with roxygen2

Exercise 2 - Create a package

From RStudio, create a New Directory (File -> New Project...)

Exercise 2 - Create a package

Exercise 2 - Create a package

  • Save functions in the R folder as .R files
  • Documentation will be automatically created in the man folder
  • You may need to delete the NAMESPACE file once to avoid warnings

Exercise 2 - use roxygen2

We start by installing roxygen2 and loading it:

install.packages("roxygen2")
library("roxygen2")

To generate documentation from our roxygen2 comments, which are denoted by the #' tags, run

roxygenise()

or press Ctrl + Shift + D.

We can now read the documentation of our function by calling

?my_function # use the name of your documented function

Exercise 2 - My solution

Testing and Debugging

Error handling

Very helpful in complex functions

subtract(2, 1)
[1] 1
subtract("2", 1)
Error in arg1 - arg2: non-numeric argument to binary operator

Check that input is correct and display custom error messages:

subtract <- function(arg1, arg2) {
  if (is.numeric(arg1) == FALSE) {
    stop("`arg1` must be numeric") 
    }
  arg1 - arg2
}
subtract("2", 1)
Error in subtract("2", 1): `arg1` must be numeric

tryCatch

Use if you anticipate an error but want function to continue.

Let’s try to generate data from a multivariate normal distribution:

mvnfast::rmvn(1, mu = c(0,0), sigma = matrix(rep(1,4), 2))
Error in mvnfast::rmvn(1, mu = c(0, 0), sigma = matrix(rep(1, 4), 2)): chol(): decomposition failed

mvnfast::rmvn is fast but fails for some problematic sigma values. In case it fails, we use MASS::mvrnorm instead:

my_rmvn <- function(n, mu, sigma) {
  tryCatch(mvnfast::rmvn(n, mu, sigma),
           error = function(e) MASS::mvrnorm(n, mu, sigma))
}
my_rmvn(n = 1, mu = c(0,0), sigma = matrix(rep(1,4), 2))
[1] -0.1206228 -0.1206228

Traceback

If something went wrong, find out where using traceback():

my_rmvn(n = 1, mu = 0, sigma = matrix(rep(1,4), 2))
Error in MASS::mvrnorm(n, mu, sigma): incompatible arguments
traceback()

Break points

Break points allow you to look inside your function’s environment. Click next to a line of code in your function to activate a break point (a red dot appears):

Now run the function…

my_rmvn(n = 1, mu = c(0,0), sigma = matrix(rep(1,4), 2))

You may need to run devtools::load_all() first.

Break points

Break points

We can now browse the function environment in the console like we normally can browse the global environment. For example we can look at sigma:

Press Stop to end the browsing. Don’t forget to deactivate the break point by clicking on the red dot in the script.

Tests

Functions should be tested before they are used.

Sometimes, interactive testing may be enough.

For example, identical() tests whether two objects are exactly equal:

identical(subtract(2,1),
          1)
[1] TRUE

Automated testing with the testthat R package is superior.


Read more on testing at r-pkgs.org.

testthat

  • Used to create unit tests for R packages that can be run automatically
  • R package, install with install.packages("testthat")
  • Set up infrastructure for a test with usethis::test_name

Create a test file for the subtract function:

usethis::use_test("subtract")

This has created a test-subtract.R file in package_name/tests/testthat/

testthat

In the test-subtract.R file, we can write tests, for example:

test_that("subtraction works", {
  expect_equal(subtract(2,1),
               1)
})

If we run this test, we should get

testthat

We can also run all of our tests at once with test_local or test_package("my_package"):

testthat::test_local()

If you have the devtools package installed, you can also use devtools::test() by pressing Test in the Build tab of Rstudio:

testthat

Useful tests include:

  • expect_equal
  • expect_true, expect_false
  • expect_error
  • expect_snapshot (for results that are difficult to describe with code)

For testing plot functions, the expect_doppelganger() function from the vdiffr package can be used:

test_that("plot function works", {
  plot_1 <- function() plot(1,2)
  vdiffr::expect_doppelganger("plot_1", plot_1)

This will create an image in the tests/testthat/_snaps/function_name directory. Upon first calling this, inspect the image to see if it is as expected. Future tests will fail if the function call in the test changes the image.

Test coverage

To check how much of your package is covered by tests, the covr package is helpful:

covr::package_coverage()

You may need to unload the package before using the covr package, e.g. with

devtools::unload()

Test coverage

The report() function lets you check test coverage line by line:

covr::report()

Exercise 3 - Test your function

  • Add tests for your function using testthat functionality
  • Add custom error messages to your function, flagging inappropriate input

Setup for testing:

install.packages(testthat)
install.packages(usethis)
usethis::use_test("function_name")

Advanced tools

apply

  • apply and related functions apply a function to elements of arrays, lists, …
  • more concise than loops, doesn’t change global environment

For example, to get the class of the first three columns of the reefs data:

reef_classes <- rep(NA,3)
for (i in 1:3) {
  reef_classes[i] <- class(palaeoverse::reefs[,i])
}
reef_classes
[1] "character" "character" "character"

Or with apply:

apply(palaeoverse::reefs[,1:3], MARGIN = 2, FUN = class)
   r_number        name   formation 
"character" "character" "character"

apply

Let’s have a look what happened there. apply(X, MARGIN, FUN) takes an array X (our reefs dataframe), and applies a function (FUN) to elements of that array, specified by MARGIN.

MARGIN = 2 indicates columns, MARGIN = 1 would be rows. So here we applied the class function to every column of reefs.

apply simplifies the output, so here it returned a vector with one element for each column.

lapply

lapply is similar to apply but for list or vector input. It returns a list for each element of the data.

data <- list(c(1, 2, 3), c(4, 5))
lapply(data, mean)
[[1]]
[1] 2

[[2]]
[1] 4.5

sapply does the same, but also tries to simplify the output; here we get a vector:

data <- list(c(1, 2, 3), c(4, 5))
sapply(data, mean)
[1] 2.0 4.5

vapply

vapply is a safer version of sapply, it requires the user to specify the anticipated class and length of the elements of the output:

data <- list(c(1, 2, 3), c(4, 5))
vapply(data, FUN = mean, FUN.VALUE = numeric(1))
[1] 2.0 4.5

Here, we specified that each element of the output should be a numeric of length 1.

Environment and scope

blank

Funders

  • European Union’s Horizon 2020 research and innovation program (MAPAS project)
    • Grant number: 947921
  • The Royal Society
    • Grant numbers: RF_ERE_210013, RGF_R1_180020, and RGF_EA_180318
  • Juan de la Cierva-formación fellowship
    • FJC2020-044836-I / MCIN /AEI / 10.13039 /501100011033
  • ETH+ Research Grant (BECCY project)
  • FAPESP Postdoctoral fellowship
    • Grant number: 2022/05697-9
  • Population Biology Program of Excellence Postdoctoral Fellowship (the University of Nebraska-Lincoln)
  • Lerner-Gray Postdoctoral Research Fellowship (American Museum of Natural History)